Nonlinear Bose-Einstein-condensate Dynamics Induced by a Harmonic 
Modulation of the s-wave Scattering Length 
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In a recent experiment, a Bose-Einstein condensate of 7 Li has been excited by a harmonic mod- 
ulation of the atomic s-wave scattering length via Feshbach resonance. By combining an analytical 
perturbative approach with extensive numerical simulations we analyze the emerging nonlinear 
dynamics of the system on the mean-field Gross-Pitaevskii level at zero temperature. Resulting 
excitation spectra are presented and prominent nonlinear features are found: mode coupling, higher 
harmonics generation and significant shifts in the frequencies of collective modes. We indicate 
how nonlinear dynamical properties could be made clearly observable in future experiments and 
compared to our results. 

PACS numbers: 03.75.Kk, 03.75.Nt, 67.85.De 



I. INTRODUCTION 

Soon after the first experimental realization of an 
atomic Bose-Einstein condensate (BEC), the understand- 
ing of the physical properties of the system has been 
tested by comparing experimental results for the frequen- 
cies of collective excitation modes [l|, with correspond- 
ing theoretical results 0-0] ■ The same paradigm remains 
of great relevance for present and future studies that fo- 
cus on a BEC in optical lattices [8J, in disordered po- 
tentials 9], for low-dimensional cold bosons [Io| or 
ultracold fermions [TJ, [l2| . The essential merit of testing 
theoretical predictions using collective modes stems from 
the possibility to measure frequencies of collective modes 
with a high accuracy of the order of 1% [l4[. For in- 
stance, a promising experimental proposal for observing 
a quantum anomaly based on a very precise measurement 
of the collective BEC mode, has been recently suggested 
in Ref. 

Usually, collective modes are induced by a modulation 
of the external trapping potential 0, 0, Il6l - [l9j . In the 
recent paper [2(j, however, an alternative way of a con- 
densate excitation has been experimentally realized. A 
broad Feshbach resonance of 7 Li [2l| allows to modu- 
late the atomic s-wave scattering length by modulating 
an external magnetic field. Based on this property, a 
harmonic modulation of the s-wave scattering length has 
been achieved 
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with a av ~ 3ao > 5 a ~ 2ao, where ao is the Bohr radius, 
yielding a time-dependent interaction among atoms. For 
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the specific set of experimental parameters, a quadrupole 
oscillation mode has been excited in this way and reso- 
nances located at the quadrupole frequency and its sec- 
ond harmonic have been observed. There are several ad- 
vantages of such experimental scheme: for instance, in 
future experiments with multispecies BEC, a single com- 
ponent could be individually excited in this way while 
the excitation of other components would occur only in- 
directly. 

Motivated by the latter experimental study, in this pa- 
per we analyze the condensate dynamics induced by the 
harmonic modulation of the s-wave scattering length in 
more detail. On the mean-field level at zero tempera- 
ture, the BEC dynamics is captured by the nonlinear 
Gross-Pitaevskii (GP) equation for the condensate wave- 
function and the time-dependent interaction leads to a 
time-dependent nonlinearity. Depending on the closeness 
of the external modulation frequency fi to one of conden- 
sate's eigenmodes, a qualitatively different dynamical be- 
havior emerges. In the non-resonant case, we have small- 
amplitude oscillations of the condensate size around the 
equilibrium widths, and we are in the regime of linear re- 
sponse. However, as fi approaches an eigenmode, we ex- 
pect a resonant behavior which is characterized by large 
amplitude oscillations. In this case it is clear that a linear 
response analysis does not provide a qualitatively good 
description of the system dynamics. 

First theoretical studies of collective modes of BEC 
0-01 bave explored dynamical properties in the linear 
regime of small amplitude oscillations. Certain non- 
linear aspects of a condensate dynamics induced by a 
trap modulation are given in Refs. [l 7H19j j . whereas two- 
component BECs are dealt with in Refs. [22|, [23|. We 
emphasize that very small oscillation amplitudes, which 
occur in the linear regime, are experimentally hard to 
observe. On the other hand, very large amplitude oscil- 
lations lead to a fragmentation of the condensate [2(| [HI . 
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Thus, the case in between is of main experimental inter- 
est and represents the main objective of our study. 

To this end, we use an approach that is complementary 
to the recent theoretical considerations 25 30] of a BEC 



time-dependent GP equation [36 



with harmonically modulated interaction. In Ref. [28 1 
the real-time dynamics of a spherically symmetric BEC 
was numerically studied and analytically explained using 
the resonant Bogoliubov-Mitropolsky method [3l|. On 
the other hand, in our approach in order to discern in- 
duced dynamical features, we look directly at the exci- 
tation spectrum obtained as a Fourier transform of the 
time-dependent condensate width. From this type of nu- 
merical analysis we find characteristic nonlinear prop- 
erties: higher harmonic generation, nonlinear mode cou- 
pling, and significant shifts in the frequencies of collective 
modes with respect to linear response results. In addi- 
tion, we work out an analytic perturbative theory with 
respect to the modulation amplitude capable of capturing 
many of the mentioned nonlinear effects obtained numer- 
ically. 

Nonlinearity-induced frequency shifts were considered 
previously in Ref. [HI for the case of bosonic collective 
modes excited by modulation of the trapping potential, 
and also in Ref. [32J for the case of a superfluid Fermi gas 
in the BCS-BEC crossover. Our analytical approach is 
based on the Poincare-Lindstedt method l3lL l33j - [35l ] . in 
the same spirit as presented in Refs. [l^, [32, |33j. How- 
ever, the harmonic modulation of the interaction strength 
introduces additional features that require a separate 
treatment. 

In Ref. [30l | it was predicted that a harmonic modu- 
lation of the scattering length leads to the creation of 
Faraday patterns, i.e. density waves, in BEC. Up to 
now, Faraday patterns have been experimentally induced 
by harmonic modulation of the transverse confinement 
strength (3(j[, which is studied analytically and numer- 
ically in Ref. (3tJ • In this paper we focus only on the 
nonlinear properties of low-lying collective modes and do 
not consider possible excitations of Faraday patterns. 

The paper is organized in the following way: in Sec. II 
we provide the necessary theoretical background - the 
mean-field GP equation and the time-dependent varia- 
tional approximation of its solution. In Sec. Ill we con- 
sider a spherically symmetric BEC with harmonically 
modulated scattering length. We analyze the condensate 
dynamics in detail and identify characteristic nonlinear 
features. In Sec. IV we focus on the experimentally most 
interesting axially symmetric BEC. Finally, in Sec. V we 
summarize the main points of our study and discuss the 
relevance of our results for future experimental setups. 



II. METHOD 

To study the dynamics of an atomic BEC at zero tem- 
perature, we use the mean-field description given by the 
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where tp{f,t) is the condensate wave- function normalized 
to unity. On the right-hand side we have a kinetic en- 
ergy term, an external axially symmetric harmonic trap 
potential V(r) — \mui 2 {p 2 + \ 2 z 2 ) with the trap aspect 
ratio A, and a nonlinear term arising from the mean-field 
interaction between the atoms. The interaction strength 
47rh m Na is proportional to the total number of atoms 
in the condensate N and to the s-wave scattering length 
a. A specific feature of the recent experiment [201 ] is the 
harmonic modulation of the s-wave scattering length de- 
scribed by Eq. (JTJ) , yielding a time-dependent interaction 
strength g = g(t). 

In order to obtain analytical insight into the conden- 
sate dynamics induced in this way, we use the Gaussian 
approximation introduced in Refs. [H, H|. To this end 
we assume that the condensate wave function with con- 
tact interaction has the same Gaussian form as in the 
non-interacting case, just with renormalized parameters. 
Thus, we use a time-dependent variational method based 
on a Gaussian ansatz, which reads for an axially symmet- 
ric trap: 
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is a time-dependent 
4> z {t) and <j>p(t) are 



2u z (t) 2 

where Af(t) = 7r _ 4-u p (t)~ 1 M2(t) _ i 
normalization, while u„(t), u z (t), 
variational parameters. They have straightforward in- 
terpretation: u p (t) and u z (t) correspond to the radial 
and to the axial condensate width, while 4> p (t) an (j) z (t) 
represent the corresponding phases. Therefore, the above 
ansatz describes dynamics of the condensate in terms of 
the time-dependent condensate widths and phases, while 
no center of mass motion is considered here. This is due 
to the fact that modulation of the interaction strength 
does not induce this type of motion. 

Following the variational approach introduced earlier 
in Ref. @, we insert ansatz ^ to the action yielding 
the GP equation, and extremize it with respect to varia- 
tional parameters. In this way we obtain a system of or- 
dinary differential equations that govern the condensate 
dynamics. Although here we consider a time-dependent 
interaction, this does not change the main steps in the 
derivation of variational equations. By extremizing the 
action, we first obtain a coupled system of differential 
equations of the first order for all variational parameters. 
The equations for the phases <fi p and <j> z can be solved 
explicitly in terms of widths u p and u z , yielding 

<t> p = mup/ (2hu p ), cj> z = mu z / (2hu z ) . (4) 

After inserting the above expressions into equations for 
the widths, we finally obtain a system of two differential 



3 



equations of second order for u p and u z 
refer to as a Gaussian approximation: 



which we will 



u p (t) + u p (t) - 
u,(t) + X 2 u z (t) ~ 
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= 0, (5) 
= 0. (6) 



In the previous set of equations and in all equations that 
follow, we express all lengths in the units of the charac- 
teristic harmonic oscillator length I — y^h/muip and time 
in units of ujp 1 . The dimensionless interaction parameter 

p{t) is given by p{t) = y/ffnNaty/l. 
Taking into account Eq. (JT|), we have: 



p(t) = p + qcosflt, 



(7) 



where p = yf2/itNa &v jl denotes the average interaction 
strength, q = yj2/nN8 a /l is a modulation amplitude, 
and f2 represents the modulation or driving frequency. 

To estimate the accuracy of the Gaussian approxima- 
tion for describing the dynamics induced by the harmonic 
modulation of the interaction strength, we compare it 
with an exact numerical solution of GP equation. The 
experiment in Ref. f20j ] was performed with the following 
values of dimensionless parameters: 



p=15, q = 10, A = 0.021. 



(8) 



In Fig. [TJ we plot the resulting time-dependent axial and 
radial condensate widths p rms (t) and z rms (t) calculated 
as root mean square values 
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P dp\^(p,z,t)\2p2, (9) 
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and compare them with numerical solutions of Eqs. ([5]) 
and ([6]). We assume that initially the condensate is in 
the ground state. In the variational description, this 
translates into initial conditions u p (0) = u p q, u p (0) — 0, 
u z (0) = u z q, Uz(0) = 0, where u p q and u z q are time- 
independent solutions of Eqs. ([5]) and ([6]), while in GP 
simulations we reach the ground state by performing 
an imaginary-time propagation [39]. For solving the 
GP equation @, we use the split-step Crank-Nicolson 
method [39j. It is evident that we have a good qualita- 
tive agreement between the two approaches. 

The main result obtained previously by using the 
Gaussian approximation is an analytical estimate for the 
frequencies of the low- lying collective modes |5j,l6|. I n this 
paper we consider excitations induced by a modulation 
of the interaction strength, and focus on the properties of 
the quadrupole and breathing mode. We assume that the 
external trap is stationary, thus preventing excitations of 
the dipole (Kohn) mode, corresponding to the center of 
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FIG. 1. (Color online) Time-dependent axial and radial 
condensate widths calculated as root mean square averages. 
Comparison of the numerical solution of time-dependent GP 
equation with a solution obtained by using the Gaussian ap- 
proximation for the actual experimental parameters in Eq. Q 
and n = 0.05. 



mass motion. By linearizing the Eqs. ([5]) and (j6]) around 
the equilibrium widths u p o and u z q, frequencies of both 
the quadrupole ujqq and the breathing mode ujbo were 
obtained: 
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For the repulsive interaction, the quadrupole mode has 
a lower frequency and is characterized by out-of-phase 
radial and axial oscillations, while in-phase oscillations 
correspond to the breathing mode. In the case of the 
experiment [13], Eq. (fTTj) yields: 



ujqo = 0.035375, uj B o = 2.00002. 



(12) 



We emphasize that, although based on the Gaussian 
ansatz, the variational approximation reproduces ex- 
actly the frequencies of collective modes not only for the 
weakly interacting BEC but also for the strongly inter- 
acting BEC in the Thomas-Fermi regime [3, [J There- 
fore, it represents a solid analytical description of BEC 
dynamics. 
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FIG. 2. (Color online) Condensate dynamics u(t) versus t within the Gaussian approximation for p = 0.4, q = 0.1 and several 
different driving frequencies Q. We plot the exact numerical solution of Eq. (|13f> . For off-resonant driving frequencies Q, we 
also show our analytical third-order perturbative result, as explained in Section III.B. 



However, due to the nonlinear form of the underlying 
GP equation, we expect nonlinearity-induced shifts in the 
frequencies of low-lying modes with respect to the values 
in Eq. (fTTj) calculated using the linear stability analysis. 
In particular, our goal is to describe collective modes 
induced by the harmonic modulation of the interaction 
strength in Eq. ([7]). In the case of a close matching of 
the driving frequency Q and one of the BEC eigenmodes 
we expect resonances, i.e. large amplitude oscillations. 
Here, the role of the nonlinear terms becomes crucial and 
nonlinear phenomena become visible, as we discuss in the 
next section. 



and a linear stability analysis yields the breathing mode 
frequency 
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Note that the above result for the breathing mode can be 
also obtained from Eq. (jlip if we set A = 1, u p0 — u z q = 
uq, and take into account Eq. (fT4")) . 



Numerical simulations 



III. SPHERICALLY SYMMETRIC BEC 

Using a simple symmetry-based reasoning, we conclude 
that a harmonic modulation of interaction strength in 
the case of a spherically symmetric BEC, i.e. A = 1, 
leads to the excitation of the breathing mode only, so 
that u p (t) = u z {t) = u(t). This fact simplifies numerical 
and analytical calculations and this is why we first con- 
sider this case before we embark to the study of a more 
complex axially symmetric BEC. 

Thus, the system of ordinary differential Eqs. (|5|) and 
© reduces to a single equation: 

The equilibrium condensate width uq satisfies 

«o-^3-4 = 0, (14) 



The main feature of the modulation induced dynamics 
is that it strongly depends on the value of the driving 
frequency J7. To illustrate this, we set p = 0.4, q — 
0.1 and solve Eq. (H~3]) for different values of f2. From 
the linear response theory, we have Uq = 1.08183, loq — 
2.06638 and we assume that the condensate is initially in 
equilibrium, i.e. u(0) — uq,u(0) — 0. Numerical results 
are plotted in Fig. [2] Large amplitude oscillations and 
beating phenomena are present for both 57 w luq and for 
Q « 2ujq. 

The phenomenology based on Eq. (| 13[) is more sys- 
tematically shown in Fig. [3] where we plot the oscillation 
amplitude defined as (w max — u m i n )/2 versus the driving 
frequency f2. A resonant behavior becomes apparent for 
both D, « loq and fl w 2luq. In the same figure we also 
show the expected positions of resonances calculated us- 
ing the linear stability analysis. Clearly, the prominent 
peaks exhibit shifts with respect to the solid vertical lines 
representing uiq and 2loq. As expected, a stronger modu- 
lation amplitude leads to a larger frequency shift, as can 
be seen from the inset. 



5 



•8 
g 




FIG. 3. (Color online) Oscillation amplitude (u max — «min)/2 
versus driving frequency Q for p = 0.4. In the inset, we zoom 
to the first peak to emphasize that the shape and value of a 
resonance occur at a driving frequency which differs from 
ujo depends on the modulation amplitude q. The solid vertical 
lines correspond to ujq and 2ojq. 



The curves presented in Fig. [3] are obtained by an 
equidistant sampling of the external driving frequency SI. 
In addition to the expected resonances close to too and 
2wo, a more thorough exploration of solutions of the vari- 
ational equation (fTB")) shows that other "resonances" are 
present such as, e. g. at SI ss ujq/2 and SI ss 2wo/3. This 
is further demonstrated in Fig. These "resonances" 
are harder to observe numerically, since it is necessary 
to perform a fine tuning of the external frequency. How- 
ever, they clearly demonstrate nonlinear BEC properties 
and an experimental observation of these phenomena is 
certainly of high interest. We note that the observed res- 
onance pattern of the form SI sa 2uq/ti (where n is an 
integer) arises also in the case of a parametrically driven 
system described by the Mathieu equation, for instance, 
in the context of the Paul trap [40( • 

To examine such excited modes directly, we look at the 
Fourier transform of the condensate width u{t). To this 
end, we numerically solve Eq. (|13[) and find the Fourier 
transform of its solution using the Mathcmatica software 
package [ill ] . An example of such an excitation spectrum 
for p — 0.4, q — 0.1, and SI = 2 is given in Fig. [5j The 
spectrum contains two prominent modes - a breathing 
mode of frequency uj (close, but not equal to uiq) and a 
mode that corresponds to the driving frequency SI, along 
with many higher-order harmonics which are of the gen- 
eral form mSl + nu), where m and n are integers. 

In Fig. [6] we juxtapose two zoomed Fourier spectra for 
two different driving frequencies for p = 0.4 and q = 0.2. 
On the left plot, we show zoomed spectrum for = 1. 
The vertical solid line corresponds to wo and we find the 
peak in the spectrum that lies almost precisely at this 
position. On the contrary, from the right plot of Fig. HI 
that corresponds almost to the resonant excitation = 2, 
we see that the prominent peak is displaced from the 
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FIG. 4. (Color online) Exact numerical solution of Eq. (| 13 j) 
for the condensate width u(t) versus t for p — 0.4, q — 0.3, 
corresponding to cjo = 2.06638. We observe large amplitude 
oscillations for SI w uio/2 in the top panel, while in the bottom 
panel the "resonant" behavior is present for Q « 2a;rj/3. 
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FIG. 5. Fourier transform of u(t) for p — 0.4, q = 0.1, 
and fi = 2. First plot presents the complete spectrum on a 
semi-log scale, while the subsequent plots focus on regions of 
interest in the spectrum. 



vertical line. This is the most clear-cut illustration of the 
shifted eigenfrequency arising due to the nonlinearity of 
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FIG. 6. Parts of the Fourier spectra for p = 0.4, q — 0.2, 
and two different driving frequencies: Q, — 1 (left) and Q — 2 
(right). Position of a linear response result uio is given by a 
vertical solid line. 



the underlying dynamical equations. Our objective is 
to develop an analytical approach capable of taking into 
account these nonlinear effects. 



B. Poincare-Lindstedt method 

In its essence, our analytical approach represents the 
standard Poincare-Lindstedt method [3l|, I33rl35l ]. Lin- 
earizing the variational equation (|13[) around the time- 
independent solution u Q for vanishing driving q = 0, we 
obtain the zeroth-order approximation for the collective 
mode oj = ujq expressed by Eq. (fl"5j) . To calculate the 
collective mode to higher orders, we explicitely introduce 
the sought-after eigenfrequency uj into the calculation by 
rescaling the time from t to s = tot, yielding the equation: 



uj 2 u(s) + u(s) - 



1 p 



u(s) 3 u(s) 4 u(s) 4 UJ 



q lC os— =0. (16) 



In the next step, we assume the following perturbative 
expansions in the modulation amplitude q: 

u(s) = u + qui(s) + q 2 u 2 (s) + q 3 u 3 (s) + . . . , (17) 

(18) 



0J = uj + q UJ\ + q 2 UJ 2 + q 3 UJ 3 + . . . , 



where we expand uj around ujq and introduce frequency 
shifts uj\, uj 2i ...for each order in the expansion in q. 
By inserting the above expansions into the Eq. (|16[) and 
collecting terms of the same order in q, we obtain a hi- 
erarchical system of linear differential equations. Up to 
the third order, we find: 

uj 2 Ui(s) +uj 2 ui(s) — -^cos 



Un UJ 



u 2 u 2 (s) + uj 2 u 2 (s) 



<>.- 



— 2ujquj\U\(s) fUi(s)cos h aui(s) , 

Uq uj 

w oW3(s) +u 2 u 3 (s) = 

— 2ujquj 2 ui{s) — 2f3ui(s) 3 + 2au\(s)u 2 (s) — ujfili(s) 

10 / x2 Ms 4 . , fis n ... . 

H — rUi(s) cos rit2 s cos 2ujoujiu 2 (s) 

UJ Un UJ 



'0 



These equations disentangle in a natural way - we solve 
the first one for iti(s) and use that solution to solve the 
second one for u 2 (s) and so on. At the n-th level of the 
perturbative expansion (n > 1) we use the initial con- 
ditions u„(0) = 0,u„(0) = 0. As is well known, the 
presence of the term cos s on the right-hand side of some 
of the previous equations would yield a solution that con- 
tains the secular term ssins. Such a secular term grows 
linearly in time, which makes it the dominant term in 
the expansion (1171) that otherwise contains only periodic 
functions in s. In order to ensure a regular behavior 
of the perturbative expansion, the respective frequency 
shifts uj\, uj 2 , . . . are determined by imposing the cancel- 
lation of secular terms. 

This analytical procedure is implemented up to the 
third order in the modulation amplitude q by using the 
software package Mathematica [41j . Although the cal- 
culation is straightforward, it easily becomes tedious for 
higher orders of perturbation theory. Note that it is nec- 
essary to perform calculation to at least the third order 
since it turns out to be the lowest-order solution where 
secular terms appear and where nontrivial frequency shift 
is obtained. We solve explicitely for ui(s), u 2 (s), and 
u 3 (s) and show an excellent agreement of our analytical 
solutions with a respective numerical solution of Eq. (fT5|) 
in Fig. [2] From the first-order solution ui(t) we read off 
only the two basic modes ujq and Q, while the second- 
order harmonics 2ujq, ujq — f2, ujq + Q and 20, appear in 
u 2 (t). In the third order of perturbation theory, higher- 
order harmonics uj — 20, 2uj — i7, 2uj + uj + 20, 3uj, and 
30 are also present. Concerning the cancellation of secu- 
lar terms, the first-order correction uj\ vanishes, leading 
to a frequency shift which is quadratic in q: 



uj = UJQ + 



P(O) 



+ ..., (19) 



where a = lOp/v,® + 6/uq and (3 = 10p/ul + 5/uq. 



12u 2Q uj3 (O 2 - uj 2 ) 2 (O 2 - 4uj 2 ) 

where the polynomial P(O) is given by 

P(fi) = O 4 [-2A0pu 5 o + 36it^(-4 + 3m 4 Wo)] 
+ O 2 [-llOOp 2 - 30pu (44 - 65u 4 ^o) 
+ 9u§(-44 + 127u 4 u; 2 - 44u^w 4 )] 
+ 5600p 2 u; 2 - 840pu o w 2 (-8 + 3u 4 w 2 ) 
+ 36wgw 2 (56 - 39m 4 u;^ + 8u 8 uj 4 ). (20) 

Mathematica notebook which implements this analytical 
calculation is available at our web site [42j]- 

C. Discussion 

The result given by Eq. (fT!)|) is the main achievement 
of our analytical analysis. It is obtained within a per- 
turbative approach up to the second order in q and it 
describes the breathing mode frequency dependence on 
O and q as a result of nonlinear effects. Due to the under- 
lying perturbative expansion, we do not expect Eq. (|19[) 
to be meaningful at the precise position of the resonances. 
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FIG. 7. (Color online) Frequency of the breathing mode ver- 
sus the driving frequency Q for p — 0.4 and q — 0.1 (upper 
plot), and p — 1 and q — 0.8 (lower plot). The dashed line 
represents fi/2 and is given to guide the eye. 



nances in the BEC with a harmonically modulated inter- 
action. A perturbative expansion to higher orders would 
probably introduce some additional poles, responsible for 
higher-order "resonant" behavior observed at ft ~ 2u!o/n, 
(n > 3). Still, the poles seem to be only an artefact of our 
approximative perturbative scheme, not present in the 
exact description. For example, a simple resummation 
performed by using the second-order perturbative result 
removes these effects, although this is only an ad-hoc ap- 
proximation. We stress that this issue concerning the 
true resonant behavior can not be settled either by rely- 
ing on a numerical calculation due to inherent numerical 
artefacts related to finite numerical precision and finite 
computational time. To resolve it, one should rely on an 
analytical consideration along the lines of Ref. [l7| or use 
some analytical tool applicable at resonances, such as the 
resonant Bogoliubov-Mitropolsky method [3l|. However, 
this is out of scope of the present paper. 

In addition to comparison of our analytical results with 
numerical solutions based on the Gaussian variational ap- 
proximation, we present a comparison with the full nu- 
merical solution of the GP equation. In order to be able 
to preform Fourier analysis with good resolution, it is 
necessary to obtain an accurate solution for long evolu- 
tion times. We do this by using the split-step method 
in combination with the semi-implicit Crank-Nicolson 
method [39| . As we refine the GP numerics by using 
finer space and time discretization parameters our nu- 
merical results become stable as shown in Fig. [3] From 
the same figure, we observe quantitatively good agree- 
ment between GP numerics and Gaussian approxima- 



However, by comparison with numerical results based on 
the variational equation, we find that Eq. (|19p represents 
a reasonable approximation even close to the resonant 
region. 

To illustrate this, in Fig. [7] we show two such compar- 
isons. In the upper panel we consider the parameter set 
p = 0.4 and q = 0.1 and observe significant frequency 
shifts only in the narrow resonant regions. We notice 
an excellent agreement of numerical values with the an- 
alytical result given by Eq. (fH)l) . In the lower panel we 
consider the parameter set p = 1 and q = 0.8 with much 
stronger modulation amplitude. In this case we observe 
significant frequency shifts for the broader range of mod- 
ulation frequencies fl. In spite of a strong modulation, we 
still see a qualitatively good agreement of numerical re- 
sults with the analytical prediction given by Eq. (TT91 . In 
principle, better agreement can be achieved using higher- 
order perturbative approximation. The dashed line on 
both figures represents fi/2, given as a guide to the eye. 
It also serves as a crude description of what we observe 
numerically in the range f2 ~ 2wo- 

The presence of two poles at 51 = ujq and — 2uq 
in Eq. (|T9")) implies the possible existence of real reso- 
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FIG. 8. (Color online) Part of the Fourier spectrum of the 
time-dependent condensate width for p = 0.4, q — 0.2, Q. = 2. 
For numerical solution of GP equation we use several dis- 
cretization schemes: GP numerics 1 (time step e = 10 -3 , 
spacing h = 4 X 10~ 2 ), GP numerics 2 (e = 5 x 10~ 4 , 
h = 2 x 10~ 2 ), GP numerics 3 (e = 5 x 10~ 5 , ft = 5x 10" 3 ). 
For comparison we also show the corresponding spectrum ob- 
tained from the Gaussian approximation (dotted-dashed line) 
and analytical result (|19[) for the position of breathing mode 
(solid vertical line). 
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tion, reflected in close values obtained for the breathing 
mode frequency. In addition, numerical values for the 
breathing mode approach closely the analytical result of 
Eq. (fit?]) , shown by a solid vertical line in Fig. [3] 

It is well known that for a corresponding two-dimen- 
sional axially-symmetric system with a constant inter- 
action and trapping frequency, the breathing mode os- 
cillations can be described by an exact linear equation 
[43|, However, in the case of a time-dependent trap- 
ping frequency, the exact equation of motion is nonlinear 
(17| . To the best of our knowledge, for a time-dependent 
interaction strength the corresponding exact equation 
does not exist in the literature, but one can reasonably 
expect that nonlinear effects will remain in such systems, 
due to the inherent time dependence of the interaction. 



and for n = 2 we get correspondingly 
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IV. AXIALLY SYMMETRIC BEC 

To obtain experimentally more relevant results, we 
now study an axially symmetric BEC. An illustration 
of the condensate dynamics is shown in Fig. [5] for p = 1, 
q = 0.2, A = 0.3. We plot numerical solutions of Eqs. (O 
and ^ obtained for the equilibrium initial conditions 
M p (0) = u p o, u p (0) = 0, u z (0) = Uzo, and u z {0) = 0. 
For the specified parameters, the equilibrium widths are 
found to be u p0 = 1.09073, u z0 = 2.40754 and from 
the linear stability analysis we find both the quadrupole 
mode frequency loqq = 0.538735 and the breathing mode 
frequency wbo = 2.00238. For a driving frequency O 
close to ujqq, we observe large amplitude oscillations in 
the axial direction. An example of excitation spectra is 
shown in Fig. [TD] Here, we have the three basic modes 
u>q, lub, and many higher-order harmonics. 



A. Poincare-Lindstedt method 

In order to extract information on the frequencies of 
the collective modes beyond the linear stability analysis, 
we apply the perturbative expansion in the modulation 
amplitude q: 

U p {t) = u p() + qu p i(t) + q 2 u p2 (t) + q 3 u p3 {t) + (21) 

u z (t) = u z0 + qu zl (t) + q 2 u z2 {t) + q 3 u z3 (t) + . . . , (22) 

By inserting these expansions in Eqs. ([5]) and © and 
by performing additional expansions in q, we obtain a 
system of linear differential equations of the general form: 

u pn (t) + mnUp n (t) + m 12 u zn (t) + f pn {t) = 0, (23) 
m 2 iu pn (t) + u zn (t) + m 22 u zn (t) + f zn {t) = 0, (24) 
where n = 1, 2, 3, . . ., ran = 4, m 12 = p/u 3 u 2 z0 , m 2 i = 
2 P/u 3 p0 u 2 z0 , and m 22 = A 2 + 3/u^ + 2p/u 2 p0 u 3 z0 . The 
functions f pn (t) and f zn (t) depend only on the solutions 
u P i(t) and u Z i(t) of lower order, i. e. i < n. For n = 1 
we have 



The linear transformation 



f P i(t) 



cos fit 

u 3 p0 u z0 1 



cosm 

Jzl(t) - - 2 , 



u pn (t) = x n (t) + y n (t), 
u zn {t) = ci x n {t) + c 2 y n (t), 



(25) 
(26) 



with coefficients 



ci 



C2 



m 2 2 - mn - \J {m 22 - mn) 2 + Am 12 m 21 
2mi2 

m 2 2 - mn + \J (TO22 — mn) 2 + 4mi2m2i 
2toi 2 



decouples the system at the n-th level and leads to: 



0, (27) 
0. (28) 



X n {t) +UJ Q0 X n (t) H ^ 

C 2 Cl 

■■ n\ i 2 n\ i Cl fp n W — fzn(t) 

y n (t) + uJ B0 y n {t) + 

ci - c 2 



Now it is clear how to proceed: we first solve Eqs. (|27|) 
and (|28| for si(i) and j/i (i) and then using Eqs. (|25|) and 
(121)1) we obtain u p i(i) and u 2 i(t). In the next step we use 
these solutions and solve for x 2 (t) and y 2 (t) and so on. 
At each level we impose the initial conditions u pn (0) = 0, 
^pn(O) = 0, u zn (0) — 0, and u zn (0) = 0. At the first 
level of perturbation theory, equations for x and y are 
decoupled, i.e. x\{t) and yi(t) are normal modes: x\{t) 
describes quadrupole oscillations, while yi{t) describes 
breathing oscillations. However, at the second order of 
perturbation theory yi(t) enters the equation for x 2 (t) 
and also Xi(t) appears in equation for y 2 (t), i.e. we have 
a nonlinear mode coupling. 

The explicit calculation is performed up to the sec- 
ond order by using the software package Mathematica 
[ill ]. We obtain a good agreement of analytical results 
obtained in the second order of our perturbation the- 
ory and numerical results, as can be seen in Fig. |9]for a 
moderate value of a modulation amplitude q. The first 
secular terms appear at the level n — 3. The expressions 
are cumbersome, but the relevant behavior is obtained 
from the terms: 



u pO u zO 



%3(t) + UQ Q X S (t) + Cq COs(uiQ t) 



0. 



(29) 
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FIG. 9. (Color online) Condensate dynamics within the Gaussian approximation for p = 1, q = 0.2, A = 0.3 and two different 
driving frequencies Q = 0.4 (left plot) and = 1 (right plot). We plot exact numerical solution of Eqs. ([5]) and © together 
with the analytical second-order perturbative result, as explained in Section IV. A. 
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for p = 1, q — 0.2, A = 0.3, and £2 = 0.4. Upper plot gives 
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the spectrum together with positions of prominent peaks. 



3=3(0 



Cc 



t sin(wQoi) 



(30) 



The last term can be absorbed into the first-order solu- 
tion 

Cr" 2 



u p (t) = A Q cos(ujQ t) 

PS Aq COs(t(u}Q 



2^qo 
Awqo)) 



t sin(wQoi) 



(31) 



and can be interpreted as a frequency shift of the 
quadrupole mode to the order q 2 : 



Cc 



2ujq Aq 



(32) 



The coefficients Aq and Cq are calculated using the 
Mathematica code available at our site [42|, but their 
explicit form is too long to be presented here. Along the 
same lines we calculate the frequency shift of the breath- 
ing mode. 



B. Discussion 

The main results of our calculation are shown in 
Figs. QT] and [12] In Fig. [TT] we plot the analytically 
obtained frequency of the quadrupole mode versus the 
driving f2, using the second order perturbation theory to- 
gether with the corresponding numerical result based on 
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versus driving frequency fl for p = 1, q = 0.2, and A = 0.3. 
We plot linear response result ujbo, second-order analytical 
result 0JB,(a) arm numerical values u>B,(n)- 



the Fourier analysis of solutions of Eqs. §5§ and ([6]). An 
analogous plot for the frequency of the breathing mode 
is given in Fig. [12] Our analytical perturbative result for 
the shifted quadrupole mode frequency contains poles at 
cjqo, 2ujq , ujbo - ujqq, ujqo + ujbo and ujbo- Similarly, 
for the shifted frequency of the breathing mode poles in 
the perturbative solution are found at ujqo, ujbo, Zujbo, 
ujbo — ujqo and ujqo + ujbo- In both figures we see ex- 
cellent agreement of the pcrturbatively obtained results 
with the exact numerics. 

In the experiment (20| . excitations of a highly elon- 
gated and strongly repulsive BEC were considered with 
the system parameters in Eq. ||5J). For that case, ac- 
cording to Eq. (TT2"|) we get ujqq <C ujboi an d the driv- 
ing frequency was chosen in the range (0,3u>qq). Good 
agreement of real-time dynamics obtained from the varia- 



tional approximation with the exact solution of the time- 
dependent GP simulation occurs even for long propaga- 
tion times according to Fig. [IJ which implies a good ac- 
curacy of the Gaussian approximation for calculating the 
frequencies of the excited modes. From the real-time dy- 
namics plotted in Fig. [T]we observe the excitation of the 
slow quadrupole mode as an out-of phase oscillation in 
the axial and in the radial direction. In addition, in the 
radial direction we observe fast breathing mode oscilla- 
tions. This is typical for highly elongated condensates 
45] and our analysis for the experimental parameters 
shows a strong excitation of the quadrupole mode, but 
also a significant excitation of breathing mode in the ra- 
dial direction. Due to the large modulation amplitude 
q, many higher order harmonics are excited, and, most 
importantly, we find frequency shifts of the quadrupole 
mode of about 10% in Fig [13J From the same figure 
we notice that, due to the chosen frequency range for 
fi, only resonances located at ujq and 2ujq are observed. 
The presence of nonlinear effects is already mentioned in 
Ref. [20] • However we suggest frequency shifts calculated 
here to be taken into account for extracting the resonance 
curves from the underlying experimental data. 

To achieve more clear-cut experimental observation of 
the nonlinearity-induced frequency shifts calculated in 
this paper, we suggest a different trap geometry from the 
one used in Ref. [2(| ■ Measurements of stable BEC modes 
can be performed for about 1 s and in order to extract 
precise values of the excited frequencies in the Fourier 
analysis, several oscillation periods should be captured 
within this time interval. A higher frequency of the 
quadrupole mode, that can be realized by using a larger 
trap aspect ratio A, in combination with a higher mod- 
ulation frequency would fulfill this condition. According 
to the results presented in Ref. [20| , resonant driving may 
lead to condensate fragmentation. However, our numeri- 
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cal results indicate frequency shifts of 10 % even outside 
the resonant regions according to Figs. [7] and H31 and 
this is where experimental measurements should be per- 
formed. Although an increase in A leads to a more pro- 
nounced nonlinear mixing of quadrupole and breathing 
mode and may complicate condensate dynamics further, 
it may be possible to perform a Fourier analysis of exper- 
imental data, analogous to Ref. 46], and to compare it 
with the excitation spectra presented here. To achieve a 
complete matching of experimental data and our calcu- 
lations, it may turn out that higher-order corrections to 
Eq. that arise due to nonlinear dependence of scat- 
tering length on the external magnetic field, have to be 
taken into account. 



V. CONCLUSIONS 

Motivated by recent experimental results, we have 
studied nonlinear BEC dynamics induced by a harmon- 
ically modulated interaction at zero temperature. We 
have used a combination of an analytic perturbativc 
approach, numerical analysis based on Gaussian ap- 
proximation, and numerical simulations of a full time- 
dependent Gross-Pitaevskii equation. We have pre- 
sented numerically calculated relevant excitation spec- 
tra and found prominent nonlinear features: mode cou- 
pling, higher harmonics generation, and significant shifts 
in the frequencies of collective modes. In addition, we 
have provided an analytical perturbative framework that 
captures most of the observed phenomena. The main re- 



sults are analytic formulae describing the dependence of 
collective mode frequencies on the modulation amplitude 
and on the external driving frequency for different trap 
geometries. To extend the applicability of our analytical 
approach, a perturbative expansion to higher order has 
to be performed, or an appropriate resummation of the 
perturbative series could be applied. 

The presented results could contribute to future exper- 
imental designs that may include mixtures of cold gases 
and their dynamical response to harmonically modulated 
interactions, such as pattern formation induced by the 
modulation of different time dependence of the scatter- 
ing length. In addition, our results could contribute to 
resolving beyond-mean-field effects in the collective mode 
frequencies, as proposed in Refs. [13,13) an d f° r dipo- 
lar BEC in (49[. Nonlinearity- induced shifts of collective 
modes have to be properly taken into account to clearly 
delineate them from beyond-mean-field effects. 
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